home *** CD-ROM | disk | FTP | other *** search
/ MacTech 1 to 12 / MacTech-vol-1-12.toast / Source / MacTech® Magazine / Volume 10 - 1994 / 10.08 Aug 94 / Programmers' Challenge / factor.c next >
Encoding:
C/C++ Source or Header  |  1994-07-18  |  14.0 KB  |  454 lines  |  [TEXT/MPS ]
  1. /* Bob Boonstra
  2.  * Solution strategy
  3.  *   Factoring is a field which has been the subject of a
  4.  *   great deal of research because of the implications for
  5.  *   cryptography, especially techniques that depend on
  6.  *   the difficulty of factoring very large numbers.
  7.  *   Therefore, it is possible that some of these
  8.  *   algorithms could be applied to the challenge.
  9.  *
  10.  *   However, in the event that no mathematician
  11.  *   specializing in the field chooses to enter the
  12.  *   Challenge, this relatively simple solution takes
  13.  *   advantage of some of the simplifying conditions in the
  14.  *   problem statement:
  15.  *   1) the numbers are relatively small (64 bits, or
  16.  *      ~<20 digits)
  17.  *   2) the prime factors are even smaller (32 bits, or
  18.  *      ~<10 digits)
  19.  *
  20.  *   This solution depends on no precomputed information.
  21.  *   It is based on Fermat's algorithm, described in Knuth
  22.  *   Vol II, which is especially well suited to the problem
  23.  *   because it is most efficient when the two p,
  24.  *
  25.  *   Fermat's algorithm requires ~(p-1)sqrt(n) iterations,
  26.  *   where n=u*v and u~=p*sqrt(n), v~=sqrt(n)/p.
  27.  *   Other algorithms require half as many iterations, but
  28.  *   require more calculation per iteration.
  29.  *
  30.  *   Fermat's algorithm works as follows:
  31.  *   1) Let n - u*v, u and v odd primes.
  32.  *   2) Set a = (u+v)/2 and b = (u-v)/2.
  33.  *   3) Then n = uv = a**2 - b**2
  34.  *   4) Initialize a = trunc(sqrt(n)), b=0, r=a**2-b**2-n
  35.  *   5) Iterate looking for r==0, with an inner loop that
  36.  *      keeps a=(u+v)/2 constant and increases b=(u-v)/2
  37.  *      by 1 each iteration until r becomes negative.  When
  38.  *      this happens, the halfsum a is increased by 1, and
  39.  *      the difference loop is repeated.
  40.  *
  41.  *   The algorithm in Knuth uses auxiliary variables x,y
  42.  *   for efficiency, where:
  43.  *   x = 2*a+1 and y = 2*b+1
  44.  *   This works fine in most cases, but causes overflow
  45.  *   of a longword when x,y are the full 32-bits in size.
  46.  *   So we have augmented the algorithm to deal with this
  47.  *   case.
  48.  *
  49.  *   This solution also uses an efficient integer sqrt
  50.  *   algorithm due to Ron Hunsinger, and extends that
  51.  *   algorithm to 64 bits.
  52.  */
  53.  
  54. #pragma options(assign_registers,honor_register)
  55.  
  56. #define ulong unsigned long
  57. #define ushort unsigned short
  58.  
  59. #define kLo16Bits 0xFFFF
  60. #define kHiBit 0x80000000UL
  61. #define kLo2Bits 3
  62. #define kLo1Bit 1
  63.  
  64. /*
  65.  * Macros RightShift2 and RightShift1 shift a 64-bit value
  66.  * right by 2 and 1 bits, respectively.
  67.  */
  68. #define RightShift32By2(xL,xH)                            \
  69. {                                                         \
  70.   xL >>= 2;                                               \
  71.   xL |= (xH & kLo2Bits)<<30;                              \
  72.   xH >>= 2;                                               \
  73. }
  74.  
  75. #define RightShift32By1(xL,xH)                            \
  76. {                                                         \
  77.   xL >>= 1;                                               \
  78.   xL |= (xH & kLo1Bit)<<31;                               \
  79.   xH >>= 1;                                               \
  80. }
  81.  
  82. /*
  83.  * Macros Add32To64 (Sub32From64) add (subtract) a 32-bit
  84.  * value to (from) a 64-bit value.
  85.  */
  86. #define Add32To64(rL,rH, a)                               \
  87.   temp = rL;                                              \
  88.   if ((rL += a) < temp) ++rH;
  89.  
  90. #define Add2NPlus1To64(lowR,highR,a)                      \
  91.   Add32To64(lowR,highR,a);                                \
  92.   Add32To64(lowR,highR,a);                                \
  93.   Add32To64(lowR,highR,1);
  94.  
  95. #define Sub32From64(rL,rH, s)                             \
  96.   temp = rL;                                              \
  97.   if ((rL -= s) > temp) --rH;
  98.  
  99. #define Sub2NPlus1From64(lowR,highR,s)                    \
  100.   Sub32From64(lowR,highR,s);                              \
  101.   Sub32From64(lowR,highR,s);                              \
  102.   Sub32From64(lowR,highR,1);
  103.  
  104. /*
  105.  * Macros Add64 (Sub64) add (subtract) two 64-bit values.
  106.  */
  107. #define Add64(qL,qH, eL,eH)                               \
  108.   Add32To64(qL,qH,eL);                                    \
  109.   qH += eH;
  110.  
  111. #define Sub64(qL,qH, eL,eH)                               \
  112.   Sub32From64(qL,qH, eL);                                 \
  113.   qH -= eH;
  114.  
  115. /*
  116.  * Macro Square64 multiplies a 32-bit value by itself to
  117.  * produce the square as a 64-bit value.  For this solution,
  118.  * we only need to execute this macro expansion once.
  119.  */
  120. #define Square64(a,rL,rH,temp)                            \
  121. {                                                         \
  122.   ulong lohi,lolo,hihi;                                   \
  123.   ushort aHi,aLo;                                         \
  124.                                                           \
  125.   aHi = a>>16;                                            \
  126.   aLo = a;                                                \
  127.                                                           \
  128.   rL = (lolo = (ulong)aLo*aLo)&kLo16Bits;                 \
  129.   lohi = (ulong)aLo*aHi;                                  \
  130.                                                           \
  131.   temp = ((lohi&kLo16Bits)<<1) + (lolo>>16);              \
  132.   rL |= temp<<16;                                         \
  133.                                                           \
  134.   temp>>=16;                                              \
  135.   temp += ((hihi = (ulong)aHi*aHi)&kLo16Bits) +           \
  136.                              (lohi>>(16-1));              \
  137.   rH = temp&kLo16Bits;                                    \
  138.                                                           \
  139.   temp>>=16;                                              \
  140.   temp += hihi>>16;                                       \
  141.   rH |= temp<<16;                                         \
  142. }
  143.  
  144. /*
  145.  * Macros LessEqualZero64 and EqualZero64 determine if
  146.  * 64-bit (signed) values are <= 0 or == 0, respectively.
  147.  */
  148. #define LessEqualZero64(vL,vH)                            \
  149.     ( (0>(long)vH) || ((0==vH) && (0==vL)) )
  150.  
  151. #define EqualZero64(vL,vH)                                \
  152.      ((0==vL) && (0==vH))
  153.  
  154. /*
  155.  * Macro LessEqual64 determines if one 64-bit quantity is
  156.  * less than or equal to another.
  157.  */
  158. #define LessEqual64(uL,uH, vL,vH)                         \
  159.     ( (uH< vH) || ((uH==vH) && (uL<=vL)))
  160.  
  161. /*
  162.  * Function prototypes.
  163.  */
  164. ulong sqrt64 (ulong nLo,ulong nHi);
  165. void Factor64(ulong lowHalf,ulong highHalf,
  166.               ulong *prime1Ptr,ulong *prime2Ptr);
  167.  
  168. /*
  169.  * The solution ...
  170.  */
  171. /*********************************************************
  172.  *                                                       *
  173.  * Factor64                                              *
  174.  *                                                       *
  175.  *********************************************************/
  176. void Factor64(lowHalf,highHalf,prime1Ptr,prime2Ptr)
  177. unsigned long lowHalf,highHalf;
  178. unsigned long *prime1Ptr,*prime2Ptr;
  179. {
  180. register ulong x,y,lowR,highR;
  181. register ulong temp;
  182. ulong sqrtN;
  183.  
  184. /*
  185.  * Fermat's algorithm (Knuth)
  186.  *
  187.  * Assume n=u*v, u<v, n odd, u,v odd
  188.  * Let a=(u+v)/2  b=(u-v)/2  n=a**2-b**2  0<=y<x<=n
  189.  * Search for a,b that satisfy x**2-y**2-n==0
  190.  *
  191.  * NOTE:  u,v given as being < 2**32 (fit in one word).
  192.  * Therefore a,b also are < 2**32 (and fit in one word).
  193.  *
  194.  * C1: Set x=2*floor(srt(n))+1,
  195.  *         y=1,
  196.  *         r=floor(sqrt(n))**2-n
  197.  *     x corresponds to 2a+1, y to 2b+1, r to a**2-b**2-n
  198.  *
  199.  * C2: if r<=0 goto C4
  200.  *       (algorithm modified to keep r positive)
  201.  * C3: r=r-y, y=y+2,
  202.  *     goto C2
  203.  * C4: if (r==0) terminate with n = p*q,
  204.  *         p=(x-y)/2, q=(x+y-2)/2
  205.  * C5: r=r+x, x=x+2,
  206.  *     goto C3
  207.  *
  208.  * This solution modifies the algorithm to:
  209.  * (1) reorder arithmetic on r to keep it positive
  210.  * (2) extend r to 64 bits when necessary
  211.  * (3) handle the trivial case where one of the primes is 2.
  212.  */
  213. /*
  214.  * Handle the trivial case with an even prime factor.
  215.  */
  216.   if (0 == (lowHalf&1)) {
  217.     *prime1Ptr = 2;
  218.     *prime2Ptr = (lowHalf>>1) | ((highHalf&kLo1Bit)<<31);
  219.     return;
  220.   }
  221. /*
  222.  * Compute truncated square root of input 64-bit number.
  223.  */
  224.   sqrtN = temp = sqrt64(lowHalf,highHalf);
  225.   Square64(temp,lowR,highR,y);
  226. /*
  227.  * Initialize r to s*s - n,
  228.  *   but calculate n-s*s to keep r positive, and fix
  229.  *   later when it is time to add x to r by calculating
  230.  *   r = x - (n-s*s).
  231.  */
  232.   Sub64(lowHalf,highHalf, lowR,highR);
  233. /*
  234.  * Handle perfect square case.
  235.  */
  236.   if ((0==highHalf) && (0==lowHalf)) {
  237.     *prime1Ptr = *prime2Ptr = sqrtN;
  238.     return;
  239.   }
  240.   y = 1;
  241.   highR = 0;
  242. /*
  243.  * Separate out the overflow case where x=2a+1 does not
  244.  * fit into a long word.
  245.  */
  246.   if ((temp=sqrtN) >= kHiBit-1) goto doLargeX;
  247. /*
  248.  * If sqrt(n) < 0x80000000, then 2*sqrt(n)+1 fits in one
  249.  * long word.  Also, n-trunc(sqrt(n))**2 < 2*trunc(sqrt(n))
  250.  * also fits in a long word.
  251.  */
  252.   x = 1+2*temp;
  253.   lowR = x - lowHalf;
  254.   x += 2;
  255.   do {
  256. C2:
  257.     if (lowR<=y) break;
  258.     lowR -= y;
  259.     y += 2;
  260.   } while (true); /* exit when r<=y */
  261. C4:
  262.   if (y==lowR) {
  263.     *prime1Ptr = (x-y-2)>>1;
  264.     *prime2Ptr = (x+y)>>1;
  265.     return;
  266.   }
  267.   lowR += (x-y);
  268.   y += 2;
  269. /*
  270.  * Fall through to modified algorithm if x overflows a
  271.  * long word.
  272.  */
  273.   if ((x += 2) < (ulong)0xFFFFFFFF-2) goto C2;
  274. /*
  275.  * Adjust x and y to guarantee they will not overflow.
  276.  * This requires some extra arithmetic to add 2*a+1 and
  277.  * subtract 2*b+1, but that is preferable to using two
  278.  * longs to represent each of x and y.
  279.  */
  280.   x>>=1;
  281.   y>>=1;
  282.   goto C3L;
  283.  
  284. doLargeX:
  285. /*
  286.  * x=2*a+1 no longer fits in 32 bits, so we sacrifice a
  287.  * little loop efficiency and let x=a.
  288.  * Likewise, we let y=b instead of 2*b+1.
  289.  */
  290.   lowR = x = temp;
  291.   Add32To64(lowR,highR,x);
  292.   Sub64(lowR,highR, lowHalf,highHalf);
  293.   ++x;
  294.   do {
  295.     if ( LessEqualZero64(lowR,highR) ) break;
  296. C3L:
  297.     Sub2NPlus1From64(lowR,highR,y);
  298.     ++y;
  299.   } while (true); /* exit when lowR,highR<=0 */
  300. C4L:
  301.   if ( EqualZero64(lowR,highR)) {
  302.     *prime1Ptr = x-y;
  303.     *prime2Ptr = x+y;
  304.     return;
  305.   }
  306.   Add2NPlus1To64(lowR,highR,x);
  307.   ++x;
  308.   goto C3L;
  309. }
  310.  
  311. /*********************************************************
  312.  *                                                       *
  313.  * sqrt64                                                *
  314.  *                                                       *
  315.  *********************************************************/
  316. /*
  317.  * sqrt_max4pow is the largest power of 4 that can be
  318.  * represented in an unsigned long.
  319.  */
  320. #define sqrt_max4pow (1UL << 30)
  321. /*
  322.  * undef sqrt_truncate if rounded sqrts are desired; for
  323.  * the factoring problem we want truncated sqrts.
  324.  */
  325. #define sqrt_truncate
  326.  
  327. /*
  328.  * sqrt64 is based on code posted by Ron Hunsinger
  329.  * to comp.sys.mac.programmer.  Modified to handle 64-bit
  330.  * values.
  331.  */
  332. ulong sqrt64 (ulong lowN, ulong highN) {
  333. /*
  334.  * Compute the integer square root of the integer argument
  335.  * n.  Method is to divide n by x computing the quotient x
  336.  * and remainder r.  Notice that the divisor x is changing
  337.  * as the quotient x changes.
  338.  *
  339.  * Instead of shifting the dividend/remainder left, we
  340.  * shift the quotient/divisor right.  The binary point
  341.  * starts at the extreme left, and shifts two bits at a
  342.  * time to the extreme right.
  343.  *
  344.  * The residue contains n-x**2.
  345.  * Since (x + 1/2)**2 == x**2 + x + 1/4,
  346.  *   n - (x + 1/2)**2 == (n - x**2) - (x + 1/4)
  347.  * Thus, we can increase x by 1/2 if we decrease (n-x**2)
  348.  * by (x+1/4)
  349.  */
  350.   register ulong lowResidue,highResidue; /* n - x**2 */
  351.   register ulong lowRoot,highRoot;       /* x + 1/4 */
  352.   register ulong half;                   /* 1/2     */
  353.   ulong highhalf,lowhalf,temp;
  354.  
  355.   lowResidue = lowN;
  356.   if (0 != (highResidue = highN)) {
  357. /*
  358.  * This code extends the original algorithm from 32 bits to
  359.  * 64 bits.  It parallels the 32-bit code; see below for
  360.  * comments.
  361.  */
  362.     highRoot = sqrt_max4pow; lowRoot = 0;
  363.     while (highRoot>highResidue)
  364.       RightShift32By2(lowRoot,highRoot);
  365.     Sub64(lowResidue,highResidue, lowRoot,highRoot);
  366. /*
  367.  * The binary point for half is now in the high order of
  368.  * two 32-bit words representing the 64-bit value.
  369.  */
  370.     lowhalf = lowRoot; highhalf = highRoot;
  371.     RightShift32By2(lowhalf,highhalf);
  372.     Add64(lowRoot,highRoot, lowhalf,highhalf);
  373.     if (0==highhalf) goto sqrt2;
  374.     half = highhalf<<1;
  375.     do {
  376.       if (LessEqual64(lowRoot,highRoot,lowResidue,
  377.                             highResidue)) {
  378.         highResidue -= highRoot;
  379.         highRoot += half;
  380.       }
  381.       if (0 == (half>>=2)) {
  382.         half = sqrt_max4pow<<1;
  383.         goto sqrt2a;
  384.       }
  385.       highRoot -= half;
  386.       highRoot >>= 1;
  387.     } while (true);
  388. sqrt2:
  389. /*
  390.  * The binary point for half is now in the lower of the
  391.  * two 32-bit words representing the 64-bit value.
  392.  */
  393.     half = lowhalf<<1;
  394.     do {
  395.       if ((0==highResidue) && (0==highRoot)) goto sqrt3;
  396.       if (LessEqual64(lowRoot,highRoot,lowResidue,
  397.                             highResidue)) {
  398.         Sub64(lowResidue,highResidue, lowRoot,highRoot);
  399.         Add32To64(lowRoot,highRoot,half);
  400.       }
  401.       half >>= 2;
  402. sqrt2a:
  403.       Sub32From64(lowRoot,highRoot,half);
  404.       RightShift32By1(lowRoot,highRoot);
  405.     } while (half);
  406.   } else /* if (0 == highResidue) */ {
  407. #ifndef sqrt_truncate
  408.     if (lowResidue <= 12)
  409.       return (0x03FFEA94 >> (lowResidue *= 2)) & 3;
  410. #else
  411.     if (lowResidue <= 15)
  412.       return (0xFFFEAA54 >> (lowResidue *= 2)) & 3;
  413. #endif
  414.     lowRoot = sqrt_max4pow;
  415.     while (lowRoot>lowResidue) lowRoot>>=2;
  416.  
  417. /*                 Decrease (n-x**2) by (0+1/4) */
  418.     lowResidue -= lowRoot;
  419. /*                 1/4, with binary point shifted right 2 */
  420.     half = lowRoot >> 2;
  421. /*                 x=1.  (lowRoot is now (x=1)+1/4.) */
  422.     lowRoot += half;
  423. /*                 1/2, properly aligned */
  424.     half <<= 1;
  425.  
  426. /* Normal loop (there is at least one iteration remaining)*/
  427.     do {
  428. sqrt3:
  429.       if (lowRoot <= lowResidue) {
  430. /*
  431.  *        Whenever we can, decrease (n-x**2) by (x+1/4)
  432.  */
  433.         lowResidue -= lowRoot;
  434.         lowRoot += half;
  435.       }
  436. /*              Shift binary point 2 places right */
  437.       half >>= 2;
  438. /*              x{+1/2}+1/4 - 1/8 == x{+1/2}+1/8 */
  439.       lowRoot -= half;
  440. /*              2x{+1}+1/4, shifted right 2 places */
  441.       lowRoot >>= 1;
  442. /*              When 1/2 == 0, bin point is at far right */
  443.     } while (half);
  444.   }
  445. #ifndef sqrt_truncate
  446.   if (lowRoot < lowResidue) ++lowRoot;
  447. #endif
  448. /*
  449.  * Return value guaranteed to be correctly rounded
  450.  * (or truncated)
  451.  */
  452.     return lowRoot;
  453. }
  454.